Wczytanie danych

# wczytanie pakietu
library("terra")

W pierwszym kroku musimy stworzyć listę plików (rastrów), które zamierzamy wczytać. W tym celu możemy wykorzystać funkcję list.files(), która jako argument przyjmuje ścieżkę do folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy wczytać (pattern = "\\.TIF$") oraz zwrócić pełne ścieżki do plików (full.names = TRUE).

# listowanie plików z katalogu
files = list.files("dane/landsat", pattern = "\\.TIF$", full.names = TRUE)
files
## [1] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF"
## [2] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF"
## [3] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF"
## [4] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4.TIF"
## [5] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5.TIF"
## [6] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6.TIF"
## [7] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7.TIF"

Kiedy utworzyliśmy już listę plików, możemy je wczytać przy pomocy funkcji rast() z pakietu terra i następnie wyświetlić metadane.

# wczytanie danych rastrowych
landsat = rast(files)
landsat # odwołanie się do obiektu wyświetla metadane
## class       : SpatRaster 
## dimensions  : 1853, 2846, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 602535, 687915, 5742525, 5798115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
## sources     : LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF  
##               LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF  
##               LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF  
##               ... and 4 more source(s)
## names       : LC08_~SR_B1, LC08_~SR_B2, LC08_~SR_B3, LC08_~SR_B4, LC08_~SR_B5, LC08_~SR_B6, ... 
## min values  :          21,          32,        1321,        1924,        6827,        6982, ... 
## max values  :       60441,       61519,       65451,       65192,       63873,       65454, ...

Możemy również skrócić lub zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.

names(landsat) # nazwy oryginalne
## [1] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1"
## [2] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2"
## [3] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3"
## [4] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4"
## [5] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5"
## [6] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6"
## [7] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7"
names(landsat) = paste0("B", 1:7) # skrócenie nazw
names(landsat) # nowe nazwy
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
# zamiana nazw
# names(landsat) = c("Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2")

Wczytanie danych wektorowych odbywa się w analogiczny sposób za pomocą funkcji vect().

# wczytanie danych wektorowych
poly = vect("dane/powiat_sremski.gpkg")
poly
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 1  (geometries, attributes)
##  extent      : 622510.6, 659710.4, 5754795, 5784772  (xmin, xmax, ymin, ymax)
##  source      : powiat_sremski.gpkg
##  coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
##  names       : TERYT
##  type        : <chr>
##  values      :  3026

Teraz możemy przygotować prostą wizualizację na przykładzie kanału bliskiej podczerwieni (near infrared; B5) oraz poligonu.

# wizualizacja
plot(landsat[[5]]) # alternatywnie: plot(landsat[["B5"]])
plot(poly, add = TRUE)

Operacje na rastrach

Zasięg naszego obszaru analizy ograniczony jest do powiatu śremskiego, natomiast scena satelitarna ma o wiele większy zasięg. W takiej sytuacji możemy dociąć rastry, dzięki czemu ich dalsze przetwarzanie będzie szybsze, a wynik obliczeń zajmie mniej miejsca na dysku. Do docinania rastrów służy funkcja crop() i jako argumenty musimy podać raster oraz wektor.

landsat = crop(landsat, poly)
plot(landsat[[5]])
plot(poly, add = TRUE)

Obszar rastra zmniejszył się. Jednak możemy zauważyć, że poza poligonem wartości nie zostały usunięte. Wynika to z faktu, że obraz zawsze docinany jest do obwiedni (bounding box), a wartości poza obiektem/poligonem w rzeczywistości są maskowane (tj. są oznaczane jako brakujące wartości). Aby zamaskować piksele poza poligonem należy użyć funkcji mask().

landsat = mask(landsat, poly)
plot(landsat[[5]])
plot(poly, add = TRUE)

Tą operację można przeprowadzić również w jednej linii kodu używając argumentu mask = TRUE w funkcji crop().

crop(landsat, poly, mask = TRUE)

W następnym kroku możemy w prosty sposób sprawdzić statystyki opisowe naszego zbioru danych.

summary(landsat)
##        B1              B2              B3              B4       
##  Min.   : 6577   Min.   :   32   Min.   : 6844   Min.   : 7686  
##  1st Qu.: 8039   1st Qu.: 8220   1st Qu.: 9093   1st Qu.: 8649  
##  Median : 8425   Median : 8614   Median : 9609   Median : 9353  
##  Mean   : 8650   Mean   : 8921   Mean   : 9903   Mean   : 9961  
##  3rd Qu.: 9044   3rd Qu.: 9370   3rd Qu.:10396   3rd Qu.:10796  
##  Max.   :19736   Max.   :21314   Max.   :24431   Max.   :27674  
##  NA's   :48391   NA's   :48390   NA's   :48390   NA's   :48390  
##        B5              B6              B7       
##  Min.   : 7240   Min.   : 7378   Min.   : 7357  
##  1st Qu.:14669   1st Qu.:12482   1st Qu.:10002  
##  Median :17462   Median :14082   Median :11417  
##  Mean   :18004   Mean   :14673   Mean   :12519  
##  3rd Qu.:21254   3rd Qu.:16616   3rd Qu.:14170  
##  Max.   :31081   Max.   :50865   Max.   :60130  
##  NA's   :48390   NA's   :48390   NA's   :48390

Jak możemy zauważyć wartości odbicia spektralnego dla naszego zbioru danych są w zakresie od kilku do kilkunastu tysięcy dla każdego kanału. Odbicie spektralne powinno być w przedziale od 0 do 1, w związku z czym nasze dane musimy przeskalować za pomocą poniższego równania:

\[x = x \cdot 0.0000275 - 0.2\]

Dla przykładu, wartość piksela w kanale bliskiej podczerwieni wynosi 15000. Używając powyższego wzoru musimy tę wartość przemnożyć przez 0,0000275 (scale factor), a następnie odjąć 0,2 (offset). Jako wynik otrzymamy odbicie o wartości równej 0,2125. Należy pamiętać, że każdy produkt/kolekcja posiada inny wzór i konieczne jest zapoznanie się z dokumentacją.

Nie ma potrzeby stosowania tego wzoru osobno dla każdego kanału w pętli, ponieważ operacje matematyczne w pakiecie terra są domyślnie stosowane dla wszystkich kanałów.

landsat = landsat * 2.75e-05 - 0.2
summary(landsat)
##        B1              B2              B3              B4       
##  Min.   :-0.02   Min.   :-0.20   Min.   :-0.01   Min.   :0.01   
##  1st Qu.: 0.02   1st Qu.: 0.03   1st Qu.: 0.05   1st Qu.:0.04   
##  Median : 0.03   Median : 0.04   Median : 0.06   Median :0.06   
##  Mean   : 0.04   Mean   : 0.05   Mean   : 0.07   Mean   :0.07   
##  3rd Qu.: 0.05   3rd Qu.: 0.06   3rd Qu.: 0.09   3rd Qu.:0.10   
##  Max.   : 0.34   Max.   : 0.39   Max.   : 0.47   Max.   :0.56   
##  NA's   :48391   NA's   :48390   NA's   :48390   NA's   :48390  
##        B5              B6              B7       
##  Min.   :0.00    Min.   :0.00    Min.   :0.00   
##  1st Qu.:0.20    1st Qu.:0.14    1st Qu.:0.08   
##  Median :0.28    Median :0.19    Median :0.11   
##  Mean   :0.30    Mean   :0.20    Mean   :0.14   
##  3rd Qu.:0.38    3rd Qu.:0.26    3rd Qu.:0.19   
##  Max.   :0.65    Max.   :1.20    Max.   :1.45   
##  NA's   :48390   NA's   :48390   NA's   :48390

Nadal możemy zauważyć, że pewne wartości przekraczają nasz zakres od 0 do 1. Są to wartości odstające, które zazwyczaj związane są z błędnym pomiarem lub nadmierną saturacją. Można ten problem rozwiązać na dwa sposoby:

  1. Zastąpić te wartości brakiem danych (NA).
  2. Dociąć do minimalnej i maksymalnej wartości.

Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.

# sposób nr 1
landsat[landsat < 0] = NA
landsat[landsat > 1] = NA
# sposób nr 2
landsat[landsat < 0] = 0
landsat[landsat > 1] = 1

Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym przypadku zamiast funkcji plot() należy użyć funkcji plotRGB() oraz zdefiniować kolejność kanałów czerwonego, zielonego oraz niebieskiego. Oprócz tego należy podać maksymalną wartość odbicia dla kanałów (w naszym przypadku scale = 1). Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować rozciągnięcie kolorów używając argumentu stretch = "lin" lub stretch = "hist".

# plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1)
plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1, stretch = "lin")

Grupowanie

# wczytanie pakietu do grupowania danych
library("cluster")

Dane do modelowania muszą zostać przygotowane w odpowiedni sposób. Modele klasyfikacyjne najczęściej na etapie trenowania wymagają macierzy lub ramki danych (data frame). Dane rastrowe można przetworzyć do macierzy przy użyciu funkcji values().

mat = values(landsat)
nrow(mat) # wyświetla liczbę wierszy
## [1] 1238760

Za pomocą interaktywnej funkcji View() możemy sprawdzić jak wygląda nasza macierz.

View(mat)

Jak można zauważyć, mnóstwo jej wartości oznaczonych jest jako brak danych (głównie są to wartości poza obszarem analizy). Zazwyczaj modele nie obsługują NA, więc musimy je usunąć. Służy do tego dedykowana funkcja na.omit().

mat_omit = na.omit(mat)
nrow(mat_omit)
## [1] 640802

Teraz przejdziemy do najważniejszego etapu analizy, czyli do wytrenowania modelu. Użyjemy prostego modelu grupowania metodą k-średnich (k-means). Ten model wymaga jedynie, aby podać z góry liczbę grup/klastrów (argument centers). Jest to algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby analiza była powtarzalna musimy ustawić ziarno losowości – set.seed().

set.seed(1)
mdl = kmeans(mat_omit, centers = 5)

W wyniku powyższej operacji otrzymaliśmy m.in.:

  1. Obliczone średnie wartości grup dla poszczególnych kanałów (mdl$centers).
  2. Wektor ze sklasyfikowanymi wartościami macierzy (mdl$cluster).

Wyświetlmy te obiekty:

mdl$centers
##           B1         B2         B3         B4        B5        B6         B7
## 1 0.02308493 0.02740529 0.04435163 0.04147431 0.1727330 0.1309803 0.07958614
## 2 0.02044635 0.02579641 0.05547638 0.03750665 0.4421821 0.1492823 0.07336897
## 3 0.08051378 0.09522015 0.13257100 0.16603768 0.2716112 0.3437177 0.32147881
## 4 0.03431384 0.04119050 0.07134604 0.06554308 0.3509450 0.2042120 0.12686408
## 5 0.04921918 0.05830302 0.08442430 0.09825275 0.2417769 0.2551757 0.19605182
head(mdl$cluster) # wyświetla pierwsze 6 elementów
## [1] 4 5 5 1 1 1

Oznacza to, że pierwszy wiersz (reprezentujący pojedyncze oczko siatki) należy do grupy 3, drugi do grupy 2, trzeci do grupy 2, itd. Kolejnym etapem jest stworzenie mapy na podstawie otrzymanego wektora z klastrami.

Na początku musimy przygotować pusty wektor składający się z całkowitej liczby pikseli rastra. Można to sprawdzić za pomocą funkcji ncell(). W naszym przypadku jest to 1 238 760.

vec = rep(NA, ncell(landsat)) # przygotuj pusty wektor

Następnie musimy przypisać nasze grupy w wektorze w odpowiednie miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych wartości można odwołać się przez funkcję complete.cases().

# zastąp tylko te wartości, które nie są NA
vec[complete.cases(mat)] = mdl$cluster

W ostatnim kroku należy skopiować metadane obiektu landsat, ale tylko z jedną warstwą, i przypisać mu wartości wektora vec.

clustering = rast(landsat, nlyrs = 1, vals = vec)

Sprawdźmy teraz jak wyglądają utworzone grupy na mapie.

colors = rainbow(5, alpha = NULL) # wybierz 5 kolorów z wbudowanej palety `rainbow`
plot(clustering, type = "classes", col = colors)

Istotą grupowania jest stworzenie grupy składających się z podobnych elementów. Natomiast naszym zadaniem jest interpretacja, co przedstawiają utworzone grupy oraz ich nazwanie. Interpretacja jest trudnym zadaniem, a często jej wyniki są niejasne. W tym celu, niezbędna jest analiza statystyk opisowych grup oraz wspomaganie się różnymi kompozycjami (w naturalnych oraz fałszywszy kolorach). Bardzo przydatna jest również wiedza o właściwościach spektralnych obiektów.

Poniżej znajduje się przykładowy wynik takiej interpretacji.

colors = c("#086209", "#fdd327", "#d9d9d9", "#29a329", "#91632b")
category = c("lasy/woda", "pola uprawne", "odkryta gleba", "roślinność", "nieużytki")
plot(clustering, col = colors, type = "classes", levels = category)

Jeśli wynik jest satysfakcjonujący, to możemy go zapisać używając funkcji writeRaster(). Taki plik można później wczytać w R lub innym programie obsługującym dane przestrzenne (np. QGIS).

writeRaster(clustering, "clustering.tif")

Wizualizacje

Do analizy i charakteryzacji grup, zamiast statystyk opisowych, mogą zostać wykorzystane wizualizacje. Największe możliwości dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”.

ggplot2 wymaga przygotowania zbioru danych do odpowiedniej postaci. Dane muszą być przedstawione jako ramka danych w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej zmiany można dokonać w prosty sposób używając pakietu tidyr.

# install.packages(c("tidyr", "ggplot2"))
library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych

Nasz zbiór danych jest całkiem pokaźny (blisko 9 mln wartości), jednak nie ma potrzeby przedstawiania wszystkich danych na wykresie. Wymaga to więcej RAM i znacząco wydłuża czas rysowania. Prawie identyczny efekt można uzyskać wykorzystując mniejszą próbkę danych. Jako przykład zobrazujmy jedynie 10 000 wartości z każdego kanału spektralnego. Do stworzenia losowej próby służy funkcja sample(). W wyniku otrzymamy indeksy wylosowanych wierszy.

idx = sample(1:nrow(mat_omit), size = 10000)
head(idx) # wyświetl 6 pierwszy indeksów
## [1] 294762 392686 640775 538191 270373 549593

Połączmy teraz wylosowane wiersze z macierzy z odpowiednimi grupami (cbind()). Następnie macierz zamienimy na ramkę danych (as.data.frame()).

stats = cbind(mat_omit[idx, ], cluster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)
##          B1        B2        B3        B4        B5        B6        B7 cluster
## 1 0.0215950 0.0250875 0.0423575 0.0363350 0.1958350 0.1312650 0.0719200       1
## 2 0.0209625 0.0296525 0.0567400 0.0480775 0.3601475 0.1481500 0.0784650       4
## 3 0.0611400 0.0707100 0.1034350 0.1209800 0.3034425 0.3146350 0.2814150       3
## 4 0.0356200 0.0396350 0.0560250 0.0620475 0.1954225 0.2082925 0.1317050       1
## 5 0.0338875 0.0392775 0.0642200 0.0594900 0.3336650 0.2240500 0.1200450       4
## 6 0.0487650 0.0554750 0.0807475 0.0867425 0.2840550 0.2358200 0.1891525       5

Jak można zauważyć, powyższe dane mają formę szeroką (każdy kanał spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu wykorzystamy funkcję pivot_longer().

stats = pivot_longer(stats, cols = 1:7, names_to = "band", values_to = "value")

Dla formalności możemy jeszcze zmienić typ danych (klastrów i kanałów) na kategoryczny (factor). W praktyce związane jest to z uproszczeniem struktury danych (przejście ze skali ilorazowej do nominalnej).

stats$cluster = as.factor(stats$cluster)
stats$band = as.factor(stats$band)
head(stats)
## # A tibble: 6 × 3
##   cluster band   value
##   <fct>   <fct>  <dbl>
## 1 1       B1    0.0216
## 2 1       B2    0.0251
## 3 1       B3    0.0424
## 4 1       B4    0.0363
## 5 1       B5    0.196 
## 6 1       B6    0.131

Struktura danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.

ggplot(stats, aes(x = band, y = value, fill = cluster)) +
  geom_boxplot()

Zmieńmy kilka domyślnych parametrów żeby poprawić odbiór ryciny.

ggplot(stats, aes(x = band, y = value, fill = cluster)) +
  geom_boxplot(show.legend = FALSE) +
  scale_fill_manual(values = colors) +
  facet_wrap(vars(cluster)) +
  xlab("Kanał") +
  ylab("Odbicie") +
  theme_light()

Zmieniając facet_wrap(vars(cluster)) na facet_wrap(vars(band)), zamiast zestawienia kanałów w poszczególnych panelach, możemy zestawić grupy.